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ABSTRACT 

This paper presents a novel method for determining the probability that a 
supernova candidate belongs to a known supernova type (such as la, Ibc, IIL, 
etc.), using its photometric information alone. It is validated with Monte Carlo, 
and both space- and ground- based data. We examine the application of the 
method to well-sampled as well as poorly sampled supernova light curves and 
investigate to what extent the best currently available supernova models can be 
used for typing supernova candidates. Central to the method is the assumption 
that a supernova candidate belongs to a group of objects that can be modeled; 
we therefore discuss possible ways of removing anomalous or less well understood 
events from the sample. This method is particularly advantageous for analyses 
where the purity of the supernova sample is of the essence, or for those where it 
is important to know the number of the supernova candidates of a certain type 
{e.g., in supernova rate studies). 

Subject headings: supernovae: general 



Introduction 



Type la supernovae, empirically established to be standardized candles, are a staple of 
experimental cosmology. A number of future cosmological probes {e.g., DES, Pan-STARRS, 
LSST, JDEM) are planning massive sky surveys that will collect very large samples of su- 
pernovae, which will have to be classified into several known types. Supernova candidates 
are most easily classified by their spectra, but for very large ground-based surveys it may 
be impractical to obtain a spectral confirmation of the supernova type for each candidate. 
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Therefore, it is imperative to develop reliable methods of supernova classif ication based on 
photo m etric inforraation a lo ne. Photometric typ i ng ha s been described in iPoznanski et al. 



fl2002[ ). iRiess et all fl2004al ). Ijohnson and Crottd (120051 ) and ISuUivan et al 



((20051), among 



others. Most of the existing methods rely on color-color or color-magnitude diagrams for 
supernova classification. 

In this paper, we propose a novel approach to the photometric typing of supernova 
candidates. This method is based on a probability derived using a Bayesian approach, and 
is well suited for extracting the maximum amount of information out of limited data. Unlike 
methods relying on color-color or color-magnitude diagrams, which require a comparison of 
the candidate's color(s) with pre-existing tables or plots - a comparison that requires a good 
understanding of the errors and assumptions that went into the making of the literature 
data - our approach calls for a calculation of a single number and automatically takes into 
account all of the information that is available for a given candidate while incorporating the 
currently best-known models for supernova behavior. 



m 



A Bayesian me t hod i n the context of supernova light curve fitting has been used 



Barris and Tonryl (120041 ): however, it was applied specifically to type la supernovae to 



deduce redshift-independent distance moduli. In this paper, on the other hand, we describe 
a probabilistic approach to typing photometrically surveyed supernovae. Such an approach 
allows the possibility of the marginalization (integration) of the unknown, or nuisance, pa- 
rameters. In contrast to a traditional calculation, this technique is simply a calculation 
of a probability and does not involve fitting or minimization. However, it does require that 
the candidate sample be well understood; or, in other words, that each candidate be one of 
a number of hypothesized objects whose behavior can be well modeled. 

The purpose of this paper is twofold. Apart from introducing a methodology that can 
be easily applied "as is" or extended as needed, we would like to test the extent to which 
applying the best currently available supernova models helps in typing supernova candidates. 
There is no doubt that within the next several years new and improved supernova models 
will be constructed; when they are, they can be easily worked into the method. 

The paper is structured as follows. In Section [21 we derive the probability that a given 
candidate is a Type la supernova. In Section [3l we discuss how well the method works when 
applied to poorly sampled, Hubble Space Telescope GOODS data, and to well-sampled, 
ground-based SNLS data. We suggest further improvements to the method in Section 13. 5[ 
discuss its possible application to "anomalous" objects in Section 13. 6t and to fitting for 
supernova parameters, in Section 13. 7[ Conclusions are given in Section HI 
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2. The Probability 

Let us suppose that we have a sample of supernova candidates where it has been estab- 
lished that every candidate is consistent with some type of astronomical object with known 
or well-modeled photometric behavior (one approach to making sure that this is indeed the 
case is discussed in Section ESD- In our example, w e consider photome tric models for type 



Ibc IIL, IIP, Iln, and standard, or "Branch- normal" (IBranch et al.lll993l ). type la supernovae 
because they are currently best known; however, the method can be trivially extended to 
include other supernova types, as well as variable objects that are not supernovae, once reli- 
able models for such objects are available. We would like to determine the probability that 
a given candidate in the sample is a supernova of some known type T, given its measured 
light curve data. Here, we will focus on the case where T is la, but the method can be easily 
applied to other types as well. 

Here and for the remainder of the paper we assume that the redshift of the supernova 
candidate is perfectly known. We begin by assuming that the light curve of the candidate is 
measured in a single broadband filter a (the method can be trivially extended to any number 
of filter bands). The light curve measurements are represented in riepochs (observing times) 
as {Ai}, where i = [l,...,ncpochs]- We would like to know the probability that this candidate 
is a type T supernova, P(T\{Ai}). Technically, since {Ai} defines the candidate, we are 
calculating the probability that the type T hypothesis is true given a candidate, or: 

P{T\{Ai}) = P(T I candidate). (1) 

The probability P{T\{Ai}) depends on many parameters. We will express P(T\{Ai}) as a 
function of these parameters and, in the end, marginalize them to quantify P{T\{Ai}) . 

Assume now that we have a photometric model (which we will also refer to as a tem- 
plate), {aj}, for the expected light curve for a supernova of type T at a given redshift, 
observed in filter a. In general, this model depends on parameters such as the date of maxi- 
mum light to (or, equivalently, the tim e difference between th e dates of maximum light for the 



model and the data, tdifr ); stretch s ( jPerlmutter et all 1 19971 ). which parametrizes the width 



of the light curve (if T = la); the assumed absolute magnitude M in the restframe S-band; 
and inte rstellar extinction p arametrized, e.g., by the Cardelli-Clayton-Mathis parameters i?„ 



and A^ ( ICardelli et all 119981 ). We will refer to this collection of parameters as 

e = {tdis,s,M,R,,A,). (2) 
The parameters 9 uniquely define {aj} for type T: 

(3) 
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In this study, we would like to determine the probabihty P(Ia|{74j}) that the candidate 
is a type la supernova given the measurements {Ai}. Unfortunately, this probability is not 
directly calculable, but one can easily calculate P({ylj}|Ia), the probability of obtaining the 
measurement given that the candidate is a type la supernova. Bayes' theorem allows us to 
relate these two quantities. Naively, we may write 

^^^^ '^^'ErPiiA.mPiT) 

where P{{Ai}\Ia) is the probability to obtain data {Ai} for supernova type la, -P(Ia) contains 
prior information about type la supernovae, and the denominator is the normalization over 
all of the known supernova types T. This sum is over a finite number of supernova types, 
which gives a legitimate probability because we assume that each candidate must be one of 
the types summed in the denominator. Of course, one would like the denominator to include 
a model for every possible object that could mimic or be a supernova; in practice, we have to 
limit ourselves to a subset of major supernova types with well-known photometric models. 

When calculating P(Ia|{74i}) one must consider that a range of possible values is allowed 
for the stretch, extinction, and other parameters that characterize a Type la supernova. 
Therefore, we express P(Ia|{74i}) as a function of these parameters and then marginalize 
them: 

P{M{A,}) = J2PiOM{M) 

e 

Pi{A,}\e,la)Pie,la) 
Er j:eP{{A,}\e,T)Pi9,Ty 



(5) 



The prior probability P{9,T) contains all known information about the distributions of tdiflfj 
M, s, and other relevant parameters for type T supernovae. Note that P{6,T) could also 
include the relative probabilities of obtaining a type T supernova, i.e., the relative rates of 
the different supernova types. However, in this study we will assume no prior knowledge of 
the relative supernova rates, and instead assume the rates are all the same. 

Because 9 and T uniquely define {%}, P{{Ai}\9,T) can be written as: 

P{{Ai}\9,T)^P{{Ai}\{aj}). (6) 

The measured light curve, {A^}, is taken to be in units of counts/second, and can fluctuate 
from the model {%} according to Gaussian statistics: 

"^epochs 25 

P{{A,}\9,T) ^ P{{A,}\{a,} = n Aa,. (7) 
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In this expression, Ai and 5Ai are experimental measurements and errors for epoch i, {aj} 
represents the predicted hght curve for a type T supernova, and Aa^ is an increment of aj. 

We use t he m odels (represented by {aj}) for supernova la, Ibc, IIL, IIP, and Iln 



from iNugentl (120061 ). Several issues are of the essence here. First, the model for type la 
supernovae extends both into the UV (below 3460 A in the supernova rest frame) and 
into the IR (above 6600 A in the supernova rest frame) regions. Both of these regions are 
poorly c onstrained with t he currently available data. This is in fact the reason why some 



authors (iGuy et al.l l2005l . for example) choose to limit themselves to only the well-known 
part of the type la supernova spectrum. Second, the behavior of non-type la supernovae, 
especially at high redshifts, is not well known; the currently available models may or may 
not be an adequate representation of these objects. In our work, however, we are interested 
in exploring the extent to which the best possible supernova models currently available offer 
a discriminating power. A Bayesian approach that allows one to easily work in the uncer- 
tainties on any prior knowledge of supernovae appears quite natural for the situation. The 
supernova models used in this work are discussed in more detail in Appendix A. 

In order to calculate the prior P{9\T), we make the simphfying assumption that Ry and 
Ay, tdiff, M and s are all independent, allowing us to factorize the probability: 

P{e\T)=P{Uis\T)P{M\T)P{s\T)P{Ry,Ay\T)P{T). (8) 

Let us consider each of the terms in Eqn. [S] in turn, starting with the probability P{T) 
of observing a supernova of type T. Since we assume no prior knowledge of the relative 
rates of the various types, each type has an equal probability of appearing in the candidate 
sample. Therefore, P{T) is a constant: 

P(T) ^ ^. (9) 

where Nt is the number of supernova types considered. 

We also assume a flat prior for the difference in the dates of maximum light between the 
model and the data, tdm- In practice, we compare the measured and modeled light curves 
shifting their relative dates of maximum by one day. Marginalization over t^is thus means 
summing over a finite number of such shifts. As each shift has an equal probability and the 
prior must be normalized to 1, 



p(t^,s\T) = ^= ^^^'"^ . . (10) 

V uiu \ J J, J- A.max -uvam ^ ' 



where Atdifi is an increment of tdifr? is the number of shifts, and the maximum t"^^ and 
minimum t^Jg" set the limits on tdiff- We take Ai'diff to be 160. 
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There are potentially large uncertainties associated with how well the luminosity func- 
tion of a particular type supernova is known. This is accounted for in the prior P{M\T): 



(M-M)-' 
P 2SM'-' 

P{M\T) = ^= AM. 



2'k5M 



(11) 



Here, AM is an interval in M. We extract the mean magnitude s M in the restfr a me B - 
band from iNugentI (120061 ). and the standard deviations, 5M, from [Richardson et all (120021 ). 
Table [T] summarizes the values used. 



-7- 



Suj)eriio\'a 'iype 


M 


()M 


la 


-19.05 


0.30 


Ibc 


-17.27 


1.30 


IIL 


-17.77 


0.90 


IIP 


-16.64 


1.12 


Iln 


-19.05 


0.92 



Table 1: Mean magnitudes (Vega, restframe S-band) and standard deviations for different 
supernova types. 
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Note that there is an add itional dispersion in magnitudes due to lensing effects; see, for 
example, lAldering et all (120061 ) . We will ignore this effect in our study. 



Concerning the prior on the stretch s, P{s\T), we make a simplifying assumption that 
for type la supernovae the distribution for the stretch is the same for all filter bands. The 
concept of stretch itself, strictly speaking, is only defined in the restframe B- an d V- bands, 
althou gh an extension of the concept into the restframe J-band is discussed in iNobili et al 
( I2OO5I ): an implementation of this extension would introduce several extra parameters into 
culat ion. For type la's, the stretch is known to have a distribution given 



the likelihood ca 
in 



Sullivan et all (120061 ). Approximating this distribution by a Gaussian: 



P(s|Ia) 



2n6s 



■As, 



(12) 



where As is an increment in s, and we obtain s = 0.97 and 6s = 0.09. For non-type la's, the 
stretch is not defined. Practically speaking, this means that the light curves of non-type la 
candidates do not depend on the value of the stretch. We therefore assume a fiat prior such 
that: 

P{s\T) = As (for T ^ la) (13) 

where Smax and Smin can be arbitrary, provided they are consistent with the limits on the 
stretch considered for the type la case. Appendix B details why it is necessary to retain the 
stretch parameter for non-type la supernovae. 



Lastly, we parametrize the effects of the interstellar extinction as in lCardelli et all (Il998l ). 
The model supernova light curves are generated with some particular values of the Cardelli- 
Clayton-Mathis parameters and R^. Introducing the actual distributions for and 
is difficult due to lack of a generally accepted model. For example, there have been several 
indications that the Milky Way value of = 3.1 may not be generally applicable . Values in 
the range of ~2 to 3.5 have been suggested (jPatil et alll2006l : IValencic et alll2004l ). Likewise, 
there is no consensus for the distributions of A^^s, although there do e xist studies offering 
various non-analytic parametrization (for example, iHatano et all (119981 ) introduce a model 
of the extinction distribution as a function of the galaxy inclination, for both la and core- 
collapse supernovae). In our study, we compromise by considering a case of no extinction, 
as well as two cases of extinction with a moderate value of Ay = 0.4 and two different values 
of Ry, 2.1 and 3.1. All three cases (A^„ = 3) are considered equally possible. In other words, 
we take: 

P(i?.,A„|T) = i-, (14) 



where Ny are the number of possible discrete values of Ry and Ay considered. It is worth 
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noting that once more generally accepted models for Ay and appear, it will be trivial to 
introduce them into the formalism. 

We can now put everything together to calculate P{{Ai}\9,T)P{9,T) for both la 
and non-la cases; for T = la: 



Pi{Ai}\9,U)P{9M 



V - 1 V - 1 V - AM _(M~M)^ . ^ As (a-^r 

^epochs . (oj-Aj)^ 



TT (15) 



and for T 7^ la: 



e 

AM (M-Mf As 



11 1 AM (M^)" 



'^epochs ^ (aj — A^)^ 



TT -A^e (16) 



. , /27rM, 
1=1 * 



Finally, we can insert Eqns. (fT5|) and (fT6i) into ([5]) to obtain P(Ia|{Aj}). Naturally, 
if there are light curve measurements in more than one bandpass, the probability can be 
calculated for all of the available filters, each with its own photometric model dependent on 
the same parameters as in the single-filter case. For multi-filter measurements, there will be 
as many products over epochs in Eqns. ( JTSll and ( JT6II as there are passbands available. 

Note that if the redshift of a supernova is not well known, it can also be introduced as 
a nuisance parameter. Likewise, the method can be trivially extended to include any other 
parameter of interest. 

3. Demonstrating the Method 

In order to demonstrate the method, we consider the following samples: 

• Poorly sampled Hubble Space Telescope (HST) GOODS data consisting of both type 
la's and non-la's. 
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• Monte Carlo events with the same poor samphng as in the GOODS data. 

• Well-sampled ground-based data consisting of all type la's. 

In all cases, we used the simulation described in Appendix A to create template light curves 
in the filter bands considered. 



3.1. The Space-based GOODS Data 

To illustrate how the method performs on space-based data when measurement epochs 
are scarce, we use the gold ( "high confidence") and silver ( "likely but not certain" ) candidates 
from the HST GOODS sample (iGiavalisco et all 12004: [Dickinson et all l2003l : iRenzini et al 



200ll ) . The gold and silver classification is described in iRiess et all (l2004bl ) , and refers to the 
level of confidence with which the type of a supernova was determined. The sample includes 
15 gold and 5 silver la and 1 gold and 6 silver core-collapse (CC) candidates. 

In addition to the GOODS data, we use a sample collected by the Supernova Cosmology 
Project (SCP) and the High-Z collaborations in the spring-summer of 2004. The latter sam- 
ple covers the North GOODS field only and consists of 4 epochs separated by approximately 
45 days. Because the data were recorded in only two filter passbands, HST ACS F775W and 
F850LP, in this section we will restrict ourselves to only using these two bands. 

We start with the data th at have been flat-f i elded and gain-corrected by the HST 
pipeline, and use MultiDrizzle (iFruchter and Hookl |2002| ) to perform cosmic ray rejection 
and to combine dithered observations. We search for and perform simple aperture pho- 
tometry on the supernova candidates in each of the five GOODS epochs. To obtain the 
multi-epoch photometry for the GOODS North data, we combine all four epochs of the 
SCP/High-Z sample and then subtract these data from each of the North GOODS epochs 
in turn. For the GOODS South data, we combine South epochs 4 and 5, as well as 1 and 2; 
we then subtract the two combined samples separately from each of the five South GOODS 
epochs. 

We require that there be at least three data points with the signal-to- noise ratio S/N 
greater than 2, and that there be at least two data points with S/N greater than 3. This 
eliminates the "single-epoch" candidates 2003es, 2003eq, 20021g, 2003ak, 2003eu, 2003al, 
and 2002fx (la's), and 2003et and 2003er (CC's). We also eliminate one silver la candidate, 
20031v, which appears to have a residual cosmic ray contamination. Figure [1] shows the 
probabilities P(Ia|{y4j}) for the remaining gold and silver candidates. 
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Fig. 1. — P(Ia|{74j}) for the gold and silver GOODS candidates: la (red, hatched histograms) 
and CC (blue, empty histograms). 
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It is apparent that only one of the remaining type la candidates has a very low P(Ia|{y4j}). 
This candidate, 2003eb (P(Ia|{y4j}) = 2 x 10""^), appears to be more likely to be a type Iln 
than a type la. Note also that a silver la candidate, 2002fy, has a somewhat marginal 
P(Ia|{Aj}) of ~0.52, because it appears about equally likely to be a la and a Ibc. 

When discussing these results, it is important to be clear that we are not performing a 
fit but calculating a single probability, where we sum over many configurations with different 
stretches, magnitudes, etc. Nevertheless, it is useful to consider the one configuration that 
best matches the data, as this configuration will contribute the most to the final probability. 
Looking at these "best-matching" configurations is a good sanity check that the method 
is working. Figure [2] shows such best-matching configurations for two candidates, one la 
(2003bd) and one CC (2002kl). The best match for the former is a la model, while the 
best match for the latter is a Ibc model. It is apparent that the best matches describe the 
observed data well. 
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Fig. 2. — The best-matching configuration for a gold type la candidate (2003bd) (left plot) 
and for a silver type CC candidate (2002kl) (right plot). Solid lines are the templates for 
type la (left plot) and type Ibc (right plot) supernovae; filled circles with error bars are the 
data. The discontinuity of the templates in plot on the right is due to the fact that the 
available Ibc spectral template only extends to 85 days in the supernova rest frame past the 
supernova's date of explosion. 
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3.2. Monte Carlo Simulations of the GOODS sample 

It is encouraging to see that even though the data sampling for these space-based data is 
very poor, the method appears able to discriminate between different supernova types most 
of the time. In order to better understand the performance of the method on such data, we 
create Monte Carlo samples of all of the supernova types considered using the same sampling 
as that of the actual GOODS data (5 epochs separated by ~45 days). We generate Monte 
Carlo samples of type la, Ibc, IIL, Hp, and Iln supernovae, each with 500 candidates whose 
redshifts, stretches, magnitudes, and extinction parameters are selected randomly according 
to the probability distributions in Eqns. fllSp and (1161) . We simulate both ACS F850LP and 
F774W bands and impose the same signal-to-noise requirements on these simulated data as 
those we chose for the GOODS data, as described below. 

Figure [3] shows the resulting probability distributions for identifying type la supernovae 
and, as an example of the usage of the method for identifying other types, type Ibc, IIL, 
IIP, and Iln supernovae. The near-zero probabilities for type T supernovae in the P{T\{Ai}) 
plots come primarily from the events where the data populate the tails of the light curves 
and the discrimination between different types becomes particularly difficult. The uppermost 
right plot in Fig. [3] shows P(Ia|{Aj}) for the type la Monte Carlo sample as a function of 
redshift; it indicates that the lower probability events are also those with lower redshifts. 
This may simply be a consequence of the fact that our passbands of choice, F775W and 
F850LP, provide the best coverage for higher redshift candidates. 
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Fig. 3. — All but uppermost right: P(T|{Aj}) for Monte Carlo generated events. The 
hatched histograms are for (from left to right, top to bottom): T = la, Ibc, IIL, IIP, and 
Iln, and the empty histograms are for all other types T' ^ T. In all of the plots, the hatched 
and empty histograms are normalized to the area of the hatched histograms. The near-zero 
P(T\{Ai}) for type T come from the events where the data populate the tails of the light 
curves, making type discrimination difficult. The uppermost right plot shows P(Ia|{y4j}) for 
the la sample as a function of redshift. The box size is proportional to the number of events 
in the bin. 
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To select a sample of candidates of a given type, one would apply a cut, pcut, on 
P{T\{Ai}), such that P{T\{Ai}) > 

Pcut- The choice of the cut of course affects both the 
efficiency (defined as the fraction of type la candidates in a sample of la's that passed the 
cut) and the purity (defined here as the fraction of non-type la candidates in a sample of non- 
la's that passed the cut) of the sample. Using our Monte Carlo samples, we show both the 
efficiency of selecting type la's and the degree of contamination of the sample with non-type 
la's as a function of Pcut in Fig- HI (left). It is apparent that the efficiency stays roughly flat, 
while the purity of the sample increases dramatically, with higher p^ut- Figure H] (right) shows 
the efficiency and purity distributions as a function of redshift for one choice of pcnt, 0.95. 
Similarly to the tendency exhibited in Fig. [3]( uppermost right), higher redshift candidates 
are both the most efficiency identified and suffer the greatest degree of contamination. 



Monte Carlo 




Monte Carlo with p^^f = 0.95 




0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 
Redstiift 



Fig. 4. — Left plot: the efficiency for type la simulated supernovae (filled circles) and 
purity for non-type la simulated supernovae (ffiled triangles) as a function of the cut pcut on 
P{T\{Ai}) (require that P{T\{Ai}) > 

Pent)- Right plot: the efficiency for type la simulated 
supernovae (filled circles) and purity for non-type la simulated supernovae (ffiled triangles) 
as a function of redshift for pcut = 0.95. 
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It should be emphasized that demonstrating the method using Monte Carlo samples 
shows the behavior of P{T\{Ai}) in rather idealized circumstances; however, it can and 
should be used for studying the method's performance. 



3.3. The Ground-based Data 



To demonstrate the method on well-sam pled data, we use the la/Ia* candidates from 
the SNLS collaboration (lAstier et al.M2006f). H ere. "la" denotes secure la's; and "la*", 
probable la's, as defined in Howell et al.l ( 2005 ). We generate the supernova templates in 
the four MegaCam bands where the data are available for most o f the candidates: g, r, i, 
and z. Unlike the SALT fitter used by the SNLS collaboration (iGuy et al.l l2005l ). we do 
not restrict ourselves to a particular restframe wavelength range for creating the templates 
(the acceptable wavelength range considered in SALT ranges from 3460 A to 6600 A). The 
supernova templates may well be rather unreliable outside the chosen SALT range, but we 
would like to find out to what extent they still offer type discrimination. 

Figure [5] (left) shows the probabilities for the 73 SNLS la/Ia* candidates to be la's, and 
Fig. [5] (right) shows an example of the best-matching configuration with respect to a type 
la model for one of the candidates, 03D4at. 
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Fig. 5. — Left plot: The probability distribution for la (hatched histogram) and la* (empty 
histogram) SNLS candidates. The plot does not show candidates 03D3bh, 03D4ag, and 
04D3kr, which could not be identified as any known type (see Table E]). Right plot: The 
best-matching configuration for candidate 03D4at. Solid lines are the templates for a type 
la supernova; points with error bars are the data. 
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Out of the 73 SNLS la/Ia* candidates, one had a probabihty that was lower than 0.5 
and three were "undefined" {i.e., they yielded zero P{T\{Ai}) for all supernova types T 
considered) . These four candidates are summarized in Table [21 
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Candidate 


Redshift 


PM{Ai}) 


SALT P(x^|DOF) 


Comment 


OSDlgt 


0.55 


1.20e-9 


0.86 


No £/-band SALT fit 


03D3bh 


0.25 


undefined 


3.64e-17 


No z-band SALT fit 


03D4ag 


0.28 


undefined 


3.34e-27 


No z-band SALT fit 


04D3kr 


0.34 


undefined 





No z-band SALT fit 



Table 2: The four la/Ia* SNLS candidates that had low/undefined P(Ia|{Aj}). The can- 
didates with undefined probabihties are those which had zero P{T\{Ai}) for all supernova 
types T considered. Also listed are the probabilities given DOF degrees of freedom, 
P(x^|DOF), from the SALT fits. Note that for all of these candidates SALT did not fit one 
band {g or z) as the mean wavelength corresponding to the omitted band in the supernova 
restframe was outside the acceptable SALT range of [3460 A - 6600 A]. 
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It should be noted that all four of these candidates had one band that was not fit by 
SALT, as the mean wavelength corresponding to the omitted band in the supernova restframe 
was outside the acceptable SALT wavelength range. Most other candidates in Fig. O (left) 
have a near-one P(Ia|{y4j}). Only one of them, 03Dlfq, is somewhat of an outlier with a 
P(Ia|{Aj}) = 0.73. This candidate has data in only two out of the four bands, i and z. The 
four failures and the fact that for all of them there was at least one band that was outside 
the [3460 A - 6600 A] restframe range may well signal a failure of the template to provide 
an adequate description of the supernova behavior outside this range; however, the fact that 
only 4 out of 73 candidates failed to yield desired discrimination is encouraging. 

It is also interesting to compare our probabilities P{T\{Ai}) with the results of the SALT 
fits - namely, the probabilities given DOF degrees of freedom, P(x^|DOF). Table [2] shows 
that the three "undefined" candidates also had exceptionally low probabilities from the 
SALT fits. Figure [6] (left) presents the probability of the SALT fits vs. P(Ia|{Ai}). One 
feature of this plot is the spread of the x^ probabilities for the candidates for which we 
calculated P(Ia|{Aj}) to be essentially one. 
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SNLS 




0.4 0.6 O.i 
P(la|{AJ) 



SNLS 




0.4 0.6 O.f 

Redshift 



Fig. 6.— Left: The SALT fits P(x^|DOF) probability vs. P{la\{Ai}) for the 70 SNLS 
la/Ia* candidates (the figure excludes the three candidates with undefined P(Ia|{y4i}) listed 
in Table [21). The box size is proportional to the number of events in the bin. Right: the 
efficiency of selecting type la supernovae for pcut = 0.98 (we require P(Ia|{Aj}) > pcut) as a 
function of redshift. 
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Figure[6] (right) shows the efficiency of selecting type la supernovae in the SNLS sample 
for one choice of Pcut = 0.98 (where we require P(Ia|{ylj}) to be greater than Pcut), as a 
function of redshift. It is apparent that the efficiency remains fiat within errors. We are 
unable to test the purity of this selection, as the non-type la candidates from the SNLS 
collaboration are not yet publicly available; however, we will certainly be able to do so when 
the data are released. 



3.4. Systematic Effects 

In order to investigate how changing the assumptions on the priors affects our results, 
we performed a number of simple tests. For all of the tests, we varied only one prior while 
keeping the rest of them unchanged. 

• Assuming that the stretch parameter for all type la candidates is 1 increases the number 
of type la/Ia* SNLS candidates that have low (< 0.5) or undefined P(Ia|{Aj}) from 
4 to 6. It also lowers the probability for one gold la GOODS candidate, 2002hp, from 
0.83 to 0.32; and for one silver la GOODS candidate, 2002fy, from 0.52 to 0.03. 

• Assuming no extinction increases the number of SNLS candidates with low/undefined 
probabilities from 4 to 13. This change does not significantly alter the probabilities for 
the gold GOODS la candidates; however, it does lower P(Ia|{Ai}) for 2002fy to 0.15; 
and increases P(Ia|{74j}) for one silver CC candidate, 2002kl, from 0.01 to 0.80. 

• Using a fiat prior on the magnitudes does not change the number of undefined prob- 
ability SNLS candidates and moves candidate's 03Dlfq probability from 0.73 to 0.43. 
For the GOODS candidates, it lowers the probability for one silver la, 2002fy, from 
0.52 to 0.14, but does not change the CC candidates' probabilities. 

• The assumption of the fiat prior on the supernova rates is undoubtedly incorrect; 
however, it is difficult to suggest a plausible alternative since the rates for different 
supernova types (particu larly for non-type l as) are not well known, especially at high 



redshifts. Recent work (IDahlen et al.ll2004l . for example) seems to indicate that the 



rates of CC supernovae with respect to type la's are higher by approximately a factor 
of 3 at redshifts up to z ~ 1. Assuming this ratio holds at all redshifts, and modifying 
the prior accordingly, 5 SNLS candidates end up with low or undefined probabilities. 
As far as GOODS candidates are concerned, this change lowers P(Ia|{y4j}) for 2002fy 
from 0.52 to 0.27; but does not alter the distribution for the CC candidates. 
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The choice of priors in any Bayesian method should be as complete as possible, reflecting 
the general belief about the relative probabilities of various choices of parameters. For the 
data we have considered, the choice of the extinction prior appears to have the biggest effect, 
although it is clear that all of the priors can have a dramatic effect on P{Ia\{Ai}) for some 
candidates. A change in any prior appears to be particularly consequential for the candidates 
whose classification was not strong to begin with {e.g., the starting P(Ia|{Aj}) with all of 
the default priors for the most affected gold la, 2002hp, was only about 0.83; and for the 
most affected silver la, 2002fy, it was 0.52). 

3.5. Further Improvements 

While the method described in the paper appears to be quite promising, there are a 
number of improvements that are still to be explored. For instance, we have ignored the 
correlations between supernova parameters by assuming the factorizable prior in Eqn. [HI 
Provided that these correlations are well known, they can be included in a Hessian matrix 
within the Gaussians. Another improvement would be to use the Poisson distributions in 
lieu of the Gaussians in Eqn. [3, especially for low photon statistics. In general, Poisson 
distributions have the added advantage of integrating nicely to Gamma functions under 
certain circumstances. Here, the photon rates are average counts over time, so they lend 
themselves better to Gaussian approximations. However, one could just as easily use absolute 
number of counts taken for a given epoch and use Poisson errors. Using better models, both 
for the supernova behavior and for the behavior of the parameters that affect supernova 
observations {e.g., the interstellar dust), is another issue that is important for providing a 
better discrimination. The proposed method is of course only as useful as the ability of 
the models it relies on to accurately represent the actual diversity of the observed supernova 
sample; future ground- and space- based samples will yield tremendous improvements in this 
regard. 

3.6. Purging Anomalies 

The calculation of P(Ia| candidate) is highly dependent on characterizing the full range 
of objects that may be in the candidate sample. In practice, we must always consider 
the possibility that the candidate sample may be contaminated with "anomalies" that are 
merely mimicking a supernova signal - e.g., certain known variable astronomical types, 
such as AGNs, or residual "objects" resulting from poor image processing. Furthermore, 
while there are now spectral templates for many special supernova types {e.g., la 1991bg, 
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1991T, high velocity Ibc), there always exists the possibility that there are other, as of yet, 
undiscovered, supernova species. However, Eqns. [15] and [16] can be used to discard such 
anomalous supernovae with the following prescription. 

First, we define a likelihood: 

/(^,r I candidate) = P{{A,}\e,T)P{e,T) (17) 

so that 

/(T I candidate) = /( 6*, T| candidate). (18) 
e 

Then, for each candidate, we calculate the likelihoods for the various types T, /(T| candidate), 
as described above. We also generate large samples of Monte Carlo events for each type 
considered (la, Ibc, IIL, etc.), and calculate a corresponding likelihood distribution 
/(T I simulated type T candidate) for each type T. We define parameter Q such that if the 
probability of obtaining /(T| simulated type T candidate) < /(T| candidate) is smaller than 
Q, then the candidate is designated an anomaly. Effectively, Q is the confidence level that 
the candidate is not a type T supernova. Note that it is advantageous to consider each type 
T in turn, rather than the entire denominator of Eqn. [S] since this approach does not force 
one to make any assumptions about the relative numbers of supernovae of type T in the 
likelihood distributions. 

The parameter Q must be larger than the inverse of the number of simulated events for 
type T. In practice, the amount of CPU time limits the number of simulated supernovae 
that one can generate, and therefore constrains Q. Note that this method does not remove 
events whose light curves are similar to known types; Q essentially defines how different 
the light curves of an anomaly must be from the known supernova light curves for it to be 
tagged an anomaly. Large values of Q will possibly allow fewer anomalies that look similar 
to known supernova types to contaminate the "purified" sample, but will also decrease the 
statistics in the number of candidates of interest. A small Q will require that anomalies 
differ more dramatically from the known supernova types, but it increases the statistics for 
the candidates that we want to measure. Small Q values are also difficult to obtain as Q is 
constrained by the number of simulations. 

Using parameter Q to make a cut on /(T| candidate) is quite similar to using the more 
traditional analysis for candidate selection. The probability quantifies how likely it is 
that a measured distribution would have an equal or greater deviation from the model than 
the one observed. Likewise, we calculate the likelihood /(T| candidate) and find how likely 
it is to for a particular candidate's data to fluctuate to model T; if it is unlikely {i.e., the 
likelihood is less than Q) for all types, we call such a candidate an anomaly. 
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Note that this method would not effectively remove poorly sampled anomalies, as their 
probability would be relatively high for every type. The practicality of this approach for 
purging anomalies will be explored in a future paper. 

Finally, note that using a method based on this simple likelihood criterion to provide the 
actual classification of supernovae is very different from our Bayesian classification approach, 
since it does not account for the possibility that a number of supernova types might have 
the same values for 9. 



3.7. Fitting with the Maximum Likelihood 



The Bayesian Adapted Template Match (BATM) method described in lBarris and Tonry 



(120041 ) uses a maximum likelihood fit to extract supernova parameters. For completeness, 
we also describe the procedure for maximizing the likelihood for a given parameter, as well 
as for estimating the errors on the parameter. 

The formulation of l{9, T| candidate) lends itself to a maximum likelihood fit to tdiflfj M , 
s or and provided that the supernova type T is known. Indeed, we have maximized 
Z(6', T|candidate) to find the "best-matching" configurations of these variables in Section [3] 
(see Figs. |2] and [S](right plot)), where we see that it does an excellent job. However, these 
configurations were found by maximizing all the parameters at the same time. To find the 
best estimate for one parameter, we must assume the possibility that the others have some 
range of values - each with a corresponding probability. That is, we need only maximize 
one parameter at a time and marginalize the rest. In this scheme, errors are estimated by 
integrating the likelihood to obtain the upper and lower bound of the 68% confidence region. 

For instance, suppose that we want the best estimate for the magnitude, M. We first 
marginalize all the parameters in the likelihood to obtain /(M, T| candidate). Then, assuming 
that maximizing /(M, T| candidate) gives the best estimate for the fitted parameters, we 
maximize /(M, T| candidate) to obtain the best estimate for M, M^ax- With the best estimate 
calculated, we then find the uncertainties by satisfying 

Q ^ /J^— -^^-/(M,T|candidate)dM ^ /aI^.+.m^ ^|candidate)dM 

/~/(M,r|candidate)(iM ^(M, r|candidate)dAf ^ ' 

where M^ax + <7m+ and M^ax — o'm- are the upper and lower bounds of the 68% confidence 
region. This procedure would then be repeated for all the parameters for which we want 
estimates, by maximizing /(s, T| candidate), /(tdis, 7" [candidate), etc. Computationally, it 
may be necessary to replace the integrals with finite sums up to values where the terms of 
the likelihood are negligible, provided there is a single peak in the likelihood distribution. 
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Note that there are circumstances where maximizing l (M, T| candidate) does not give 
the best estimate for the fitted parameters. As an example, iBhat et aLl (119971 ) argue that a 
correction to the parameters is needed for this to be the case. This is one crucial reason why 
the performance of the likelihood always needs to be tested with Monte Carlo calculations. 



4. Conclusion 

We introduced a novel method to determine the probability that a supernova candidate 
is indeed a certain type of supernova sing photometric information alone. The probability 
is derived using a Bayesian approach. We have tested the method on both poorly sampled 
HST GOODS space-based data and well-populated SNLS ground-based light curves, with 
good results even for data where very few epochs available. While we considered primarily 
the application of this method to identifying type la supernova candidates, it can of course 
be used to identify any other known supernova types (see Fig. [3]). The results of these studies 
show the method to be promising. 

We have assumed that the candidate sample consists entirely of supernovae. The method 
naturally incorporates a number of possible hypotheses for the supernova types, allowing one 
to introduce prior knowledge of the probability distributions for various supernova parame- 
ters, as well as to marginalize them. Because it is crucial that the sample consist of candidates 
that have reliable models {i.e., well understood supernova types such as Branch-normal la, 
Ibc, etc.), we also proposed a possible approach for eliminating anomalous or less well un- 
derstood candidates from the data. As more models for different supernova types become 
available, more candidates will pass this preliminary selection to be further considered for 
typing. 

We have used the best currently available supernova models for all of the supernova 
types considered. While aware of their limitations, we note that using these models already 
yields good discriminating power; as better supernova models become available from the 
ongoing and upcoming ground- and space- based supernova surveys, classifying supernovae 
with this approach will become proportionally more reliable. 

The method described in this paper could be useful for a number of studies. An obvious 
application is selecting a high-purity sample of type la's for further analysis - e.g., creating 
a Hubble diagram and extracting cosmology. In order to do so, one would re-parametrize 
the likelihood in terms of the distance modulus fi and maximize z, Ia|candidate) to 
obtain the best estimates of the distance modulus and the redshift of the candidate, as well 
as the errors on these parameters, as described in Section 13. 7[ A pure sample of type T 
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supernovae could also be used for calculating the rates of this class of objects, provided one 
performs careful Monte Carlo studies of the method's rejection factor for type T supernovae. 
Another potentially useful application would be adapting the method for early-time light 
curve parameters in such a way as to make it a trigger for type la supernovae in large 
supernova surveys. Other possible extensions of the method would include introducing prior 
information on supernova rates for the different supernova types, as well as building a Bayes 
factor to measure the relative likelihood that a given candidate is more or less likely to be a 
la supernova than another type. 

We are grateful to the anonymous referee for many useful comments and suggestions. 
We would like to thank Reynald Pain for providing us with the SNLS type la supernova 
photometry, and Mark Sullivan for making his stretch data available to us. We are also 
grateful to Natalie Roe, Alex Kim, Kyle Barbary, Greg Aldering, Saul Perlmutter, Eric 
hinder, and Tony Spadafora for fruitful discussions. We would especially like to thank Segev 
BenZvi for his careful reading of the paper draft and his many useful suggestions. NK is 
partially supported by the Director, Office of Science, Department of Energy, under grant 
DE-AC02-05CH11231. 



A. Simulating Supernova Templates 



In order to create the template supernova light curves in the a v ailabl e filter bands, we 
do the following. We start with a spectral template from iNugentl (120061 ). The templates 
are made primarily from the publicly available supernova data assembled in the SUSPECT 
databaseo We consider five supernova t ypes: Ib c , IIL, IIP, Iln, and Branch-normal la. The 
Branch-normal la template is based on iNugentl (120021 ) . with some additional features such 
as the smo othing out of the U V part. This template is for a stretch 1 la. T he Ibc template 
is based on Levan et al. J2005h : the IIL, IIP, and Il n templates are ba sed on iGilliland et al.l 
(119991 ). with a contribution from lBaron et al.l (120041 ) for HP's and from lDi Carlo et al. 
for Iln's. 

Using the template for a given type, we compute the expected fiux of the supernova at 
a given redshift, assuming the ACDM cosmology (^a = 0.7, Qm = 0.3, and w = const = 
-1). The templates are created for 12 different values of the peak 5-band magnitudes, in 
the range of ± 3 6M, as listed in Table [TJ For type la's, each template is also generated 
for 14 different values of stretch, ranging from 0.6 to 1.3. Three types of templates are 



http: / /bruford. nhn.ou.edu/~suspect/indexl. html 
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made: one assuming no interstellar extinction, and two with extinction parametrized by the 
Cardelli-Clayton-Mathis parameters {A^,Ry) = (0.4,2.1) and (0.4,3.1). 

The simulation enables one to simulate either a ground- or a space- based observatory. 
The generated supernova flux is convolved with the atmospheric transmission (for the ground- 
based case), as well as relevant telescope, detector, and filter transmissions (for either ground- 
or space- based case). After the observations are thus generated, they are realized by an 
aperture exposure time calculator, creating light curves in each of the relevant filter bands. 
The signal-to- noise S/N ratio is computed as: 

S_ _ RsN fsig t ^^^^ 

^ \/RsN fsig t + Rsky At + RdUpt + N^^^^ Up U^xp + A^A 

where Rsn is the supernova source rate, fsig is the fraction of the source fiux in the aperture, t 
is the exposure time, Rsky is the sky rate, A is the effective seeing area, (computed taking into 
account the atmospheric seeing for the ground-based case, the diffraction radius calculated 
for a given filter central wavelength, the detector diffusion, and the detector pixel size), Rd 
is the detector dark current rate, N^ead is the detector read noise, Up is the number of pixels 
in the aperture, Uexp is the number of exposures, and A^a is the flatfielding contribution from 
the detector inter-pixel sensitivity variations: 

^Va = {{fs^g Rsn + Rsky A + Rdnp)xt A)' (A2) 
where A is the fiatfielding error, taken to be 10~^. 



B. The Stretch Parameter and Non-type la Supernovae 

The stretch parameter is not defined for non-type la supernovae. However, the cal- 
culation of P(T I candidate) necessitates that it be formally introduced for all candidates, 
regardless of type. If it was not included, the sum over the stretch in the numerator and 
denominator of Eqn. O would weight the non-type la's by a factor of Smax — Smin- 

To show how the stretch might be introduced for non-type la supernovae, we define T' 
as some generic non-type la supernova type, and 

d' = {tdiff,M,A,,R,). (Bl) 

we then calculate Y^e^iiMW ^T')P0\T') imposing Eqn.[8j 

oo 

Y,P{{A,}\eX)P{eX) = ^5^P({A,,}|^"'',.,T')P(^"'',.,T') (B2) 

e 9' 
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= J2J2 T')p{e'\T')p{s\T')p{r). 

0, -oo 

where, not specifying anything about the stretch parameter s, we allow it to take on any 
value between — oo and +00. 

Then, for type T' supernovae, we require that the likelihood, P{{Ai}\6' , s,T') remain 
unchanged for any choice of s. That is, 

P{{A,}\9',s,T') = P{{A.}\9',T'). (B3) 

Therefore, Eqn. lBSI becomes 

00 

J2 ^' T')P0'\T')P{s\T')P{T') (B4) 

01 -00 

00 

= P{{A,]\e\ T')p{e'\T')p{r) ^ p{s\t'). 

Ql -00 

If p{s\T') is the probability density such that P{s\T') = p{s\T')As, then probability theory 
requires 

00 

Ypis\T')As = 1. (B5) 

—00 

However, Eqn. [5] remains the same if we set 

p{s\T') = ; (B6) 



as well as replace the upper and lower bounds in Eqn. lBSl with Smax and Smm, respectively; and 
use the same limits of integration, Smin and Smax, for both the numerator and denominator. 
This is the reason for the choice of stretch prior in Eqn. 
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